
	/********************************************************************/
	/*                                                                  */
	/*                                                                  */
	/*              Copyright (C), 1993  SEIKO EPSON Corp.              */
	/*              ALL RIGHTS RESERVED                                 */
	/*                                                                  */
	/*                                                                  */
	/*   file name : pow.c                                              */
	/*                                                                  */
	/*   Revision history                                               */
	/*       93/08/31    first release                 K.Shinozuka      */
	/*       93/12/17    change include file           K.Shinozuka      */
	/*       97/03/17     modified	    		   T.Katahira	    */
	/*			change function declaration		    */
	/*					to the ANSI style	    */
	/*       97/03/25    change static to const static      M.Kudo      */
	/*       03/03/13    bug fixed( pow of 1.XXX erro and               */
	/*                     easy coding miss)     Iruma M.Takeishi       */
	/*                                                                  */
	/********************************************************************/

#ifdef UNIX
#include "math.h"
#include "smcvals.h"

#else
#include <math.h>
#include <smcvals.h>
#endif

#define CHECKY 800

extern double _Int(  );
extern double _Xfrac(  );



double
pow( double x, double y )

	/********************************************************************/
	/*                                                                  */
	/*  The pow function calculates x raised to the power y.            */
	/*                                                                  */
	/*  RETURN  normal: x raised to the power y                         */
	/*          error : NaN (input error)                               */
	/*                  - x is negative and y is not an integral valuie */
	/*                  - x is zero and y is less than or equal to zero */
	/*                  Inf, -Inf ( overflow )                          */
	/*                                                                  */
	/*     This function uses the external functions,                   */
	/*            exp.c, log.c , ldexp.c.                               */
	/*                                                                  */
	/********************************************************************/

/* ALGORITHM
 * compute pow(x,y) = x^y as exp( y * log(x) )
 * 1.check the sign of pow(x,y)
 *     x is negative and y is odd number --- pow(x,y) is negative
 *                                           let x = -x and sign = -1
 *     otherwise, pow(x,y) is positive and sign = 1
 * 2.|y| >= CHECKY
 *     return( sign * exp( y * log( x ) )
 * 3.let x to be 2^xexp * xfrac  ( sqrt(2)/2 <= xfrac < sqrt(2) )
 *     pow(x,y) = exp( y * xexp * ln2 ) * exp( y * log( xfrac ) )
 * 4.separate y*xexp into integral value and fraction
 *     y * xexp = i_yx + f_yx    ( 0 <= f_yx < 1 )
 * 5.breake y into integral value and fraction
 *     y = i_y + f_y    ( 0 <= f_y < 1 )
 * 6.pow(x,y) = sign * exp( i_yx * ln2 + i_y * log( xfrac ) ) 
 *                        *  exp( f_yx * ln2 + f_y * log( xfrac ) )
 *            = sign * 2^i_yx * xfrac^i_y * exp( f_yx * ln2 * f_y * log( xfrac ) ) 
 */

/*  the architecture of double floating point
 *
 *   0 1          11                  31 32                               63 bit
 *   -----------------------------------------------------------------------
 *  | |   exponent  |                    fraction                           |
 *   -----------------------------------------------------------------------
 *
 *  |               |                   |                                   |
 *  |    12 bits          20 bits       |             32 bits               |
 *  |            lower word             |            higher word            |
 *
 *         bit    0         sign bit         (  1 bit  )
 *              1 - 11      exponent part    ( 11 bits )
 *             12 - 63      fraction part    ( 52 bits )
 *
 */
 
	/* double          x;           the base                                     */
	/* double          y;           the power                                    */
	{
	 DBL_LNG        ckarg;       /* exchange the type to check the argument      */ 
	 DBL_LNG        power;       /* pow( x, y )                                  */
	 double         i_y;         /* intetger part of y                           */
	 double         f_y;         /* fraction part of y                           */
	 double         xfrac;       /* sqrt(2)/2 <= xfrax < sqrt(2)                 */
	 long           xexp;        /* x = 2^xexp * xfrac                           */
	 double         p_xfrac;     /* xfrac^i_y                                    */
	 double         i_yx;        /* integer part of y * xexp                     */
	 double         exp2;        /* exp( frac(y*xexp) * ln2 + f_y * log( xfrac ) */
	 register int   sign;        /* sign of x, 1 ( x>0 ) or -1 ( x<0 )           */
	 int            odd;         /* y is odd number odd = -1; else odd = 1       */
	 register short shorty;      /* change double y to short int y               */
	 const static double  ln2 = 0.69314718055994530942;

/* check the argument */	
	 ckarg._D = x;
	 CHECK_ARG( &ckarg );                          /* check NaN */
	                                               /* error: error code 33 domain error */   
	 if ( ckarg.st._LL == P_INF.st._LL ) {         /* x = Inf */
	    ckarg._D = y;
	    CHECK_ARG( &ckarg );                       /* check NaN */
	                                               /* error: error code 33 domain error */
	    if ( ckarg.st._LL == P_INF.st._LL ) {      /* y = Inf */
	       errno = ERANGE;                         /* 34: range error */
	       return( P_INF._D );
	    } else if ( ckarg.st._LL == N_INF.st._LL ) {  /* y = -Inf */
	       return( 0.0 );
	    } else if ( y > 0.0 ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( P_INF._D );
	    } else if ( y < 0.0 ) {
	       return( 0.0 );
	    } else {
	       return( 1.0 );
	    }
	 } else if ( ckarg.st._LL == N_INF.st._LL ) {  /* x = -Inf */
	    ckarg._D = y;
	    CHECK_ARG( &ckarg );                       /* check NaN */
	                                               /* error: error code 33 domain error */
	    if ( ckarg.st._LL == P_INF.st._LL ) {      /* y = Inf */
	       errno = EDOM;                           /* 33: domain error */
	       return( NAN._D );
	    } else if ( ckarg.st._LL == N_INF.st._LL ) {  /* y = -Inf */
	       return( 0.0 );
	    } else if ( (i_y = _Int( y, &odd )) != y ) {  /* y is not integer */
	       errno = EDOM;                           /* 33: domain error */
	       return( NAN._D );
	    } else if ( y > 0.0 ) {
	       if ( odd == 1 ) {                       /* y is even number */
	          errno = ERANGE;                      /* 34: range error */
	          return( P_INF._D );
	       } else {
	          errno = ERANGE;                      /* 34: range error */
	          return( N_INF._D );
	       }
	    } else if ( y < 0.0 ) {
	       return( 0.0 );
	    } else {
	       return( 1.0 );
	    }
	 }

	 ckarg._D = y;
	 CHECK_ARG( &ckarg );                          /* check NaN */
	 if ( ckarg.st._LL == P_INF.st._LL ) {         /* y = Inf */
	    if ( x <= -1.0 ) {
	       errno = EDOM;                           /* 33: domain error */
	       return( NAN._D );
	    } else if ( x > 1.0 ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( P_INF._D );
	    } else if ( x != 1.0 ) {
	       return( 0.0 );
	    } else {
	       return( 1.0 );
	    }
	 } else if ( ckarg.st._LL == N_INF.st._LL ){   /* y = -Inf */
	    if ( ( x < -1.0 ) || ( x > 1.0 ) ) {
	       return( 0.0 );
	    } else if ( ( x >= -1.0 ) && ( x < 0.0 ) ) {
	       errno = EDOM;                           /* 33: domain error */
	       return( NAN._D );
	    } else if ( ( x >= 0.0 ) && ( x < 1.0 ) ) {
	       errno = ERANGE;                         /* 34: range error */
	       return( P_INF._D );
	    } else {
	       return( 1.0 );
	    }
	 }

	 if ( x == 0 ) {                          /* zero */
	    if ( y <= 0 ) {                       /* y is less than or equal to zero */
	       errno = EDOM;                      /* 33: domain error */
	       return( NAN._D );
	    } else {                              /* y is positive */
	       return( 0.0 );
	    }
	 }

	 if ( x == 1.0 ) {
	    return( 1.0 );
	 } else if( x > 0 ) {                     /* x is positive */
	    sign = 1;
	 } else if ( (i_y = _Int( y, &odd )) != y ) {   /* y is not integer */
	    errno = EDOM;                         /* 33: domain error */
	    return( NAN._D );
	 } else {                                 /* y is integer */
	    x = -x;
	    sign = odd;
	 }
	 if ( y == 0 ) {
	    return( 1.0 );
	 }

	 if ( y > CHECKY || y < - CHECKY ) {
	    power._D  = exp( y * log( x ) );
	    if ( power.st._LL > P_INF.st._LL ) {  /* overflow in exp function */
	       errno = ERANGE;                    /* 34: overflow */
	       if ( sign == -1 ) {                /* negative overflow */
	          return( N_INF._D );
	       } else {
	          return( P_INF._D );
	       }
	    } else {
	       return( sign * power._D );
	    }
	 }

	 /* >>>>> add iruma m.takeishi '03.03.13 */
	 exp2 = 0.0;
	 /* <<<<< add iruma m.takeishi '03.03.13 */

	 i_yx = 0.0;
	 xfrac = _Xfrac( x, &xexp );              /* get xfrac, xexp is stored in &xexp */
	 if ( xexp != 0.0 ) {
	    i_yx = _Int( y * (double)xexp, &odd );  /* get integer */
	    exp2 = ( y * (double)xexp - i_yx ) * ln2;  // get faraction //
	    /* >>>>> delete iruma m.takeishi '03.03.13 */
	    /*
	    shorty = (short) y;
	    */
	    /* <<<<< delete iruma m.takeishi '03.03.13 */
	 }

/* calculate sin * 2^i_yx * p_xfrac * exp( exp2 ) */
	 p_xfrac = 1.0;
	 if ( xfrac != 1.0 ) {                    /* calculate 2^shorty, exp2 */
	    /* >>>>> add iruma m.takeishi '03.03.13 */
	    shorty = (short) y;
	    /* <<<<< add iruma m.takeishi '03.03.13 */

	    f_y = y - (double)shorty;
	    if ( f_y != 0.0 ) {
	       exp2 = exp2 + f_y * log( xfrac );
	    }

	    if ( shorty < 0 ) {                   /* y is negative number */
	       shorty = -shorty;
	    }
	    while ( shorty != 0 ) {               /* y is not zero */
	       if( shorty & 1 ) {                 /* the least significant bit is HIGH */
	          p_xfrac = p_xfrac * xfrac;      /* get xfrac^shorty */
	       }
	       xfrac = xfrac * xfrac;
	       shorty = ( shorty >> 1 );  
	    }
	    if ( y < 0.0 ) {
	       p_xfrac = 1 / p_xfrac;             /* xfrac^(-shorty) = 1/xfrac^shorty */
	    }
	 }

	 p_xfrac = p_xfrac * exp( exp2 );
	 power._D = ldexp( p_xfrac, (int)i_yx );  /* 2^i_yx * p_xfrac */
	 if ( power.st._LL > P_INF.st._LL ) {     // overflow in exp function //
	    errno = ERANGE;                       /* 34: overflow */
	    /* >>>>> change iruma m.takeishi '03.03.13 */
	    /*
	    if ( sign = -1 ) {                    // negative overflow //
	    */
	    if ( sign == -1 ) {                   /* negative overflow */
	    /* <<<<< change iruma m.takeishi '03.03.13 */
	       return( P_INF._D );
	    } else {
	       return( N_INF._D );
	    }
	 }
	    return( sign * power._D );
	}




double
_Int( x, odd )

	/************************************************************************/
	/*                                                                      */
	/*  The _Int function computes the largest integral value not greater   */
	/*  than x.  If x is less than 1, odd code is 0. If return value is     */
	/*  odd, the code is -1, otherwise the code is 1. The code is stored in */
	/*  the object pointed by odd.                                          */
	/*  The argument x is positive number, not zero, Inf.                   */
	/*                                                                      */
	/*  RETURN    the largest integral value not greater than  x            */ 
	/*                                                                      */
	/************************************************************************/

/* ALGORITHM
 * 1.separate x to fraction part and exponent part
 *     x = 2^xexp * xfrac
 *       xfrac: [1,2)
 * 2.if xexp is more then 51, x doesn't have fractional bits
 * 3.if xexp is less than 0, x doesn't have integral bits
 * 4.else, x has both fractional and integral bits
 * 5.clear fractional bits
 * 6.check the lowest bit of integral part, and set the odd code
 */

	double         x;           /* positive and not Inf number         */
	int            *odd;        /* the pointer, the odd code is stored */
	{
	 DBL_LNG       xfrac;       /* fraction part of x                  */
	 DBL_LNG       ipart;       /* integral number                     */
	 register long xexp;        /* exponent part of x                  */

	 xfrac._D = x;                
	 xexp = (( xfrac.st._LL & MASK_EXP ) >> L_FRAC );    /* get exponent part  */
	                                          /* MASK_EXP = 0x7ff00000 */
	 xexp = ( xexp - ((unsigned long)NUM1_L >> L_FRAC) );    /* get off offset */
	                                          /* NUM1_L = 0x3ff00000 */
	 xfrac.st._LL = ( xfrac.st._LL &  MASK_FRAC );    /* get fraction part  */
	                                          /* MASK_FRAC = 0x000fffff */
	 if ( xexp >= N_FRAC ) {                  /* no fractional bits  */
	    ipart._D = x;                         /* set integer */
	    if ( xexp == N_FRAC ) {               /* the least significant bit is 2^0 */
	       if ( ( xfrac.st._LH & 1 ) == 0 ) { /* the least significant bit is LOW */
	          *odd = 1;
	       } else {
	          *odd = -1;
	       }
	    } else {
	       *odd = 1;
	    }
	 } else if ( xexp < 0 ) {                 /* x < 1 : all fractional bits */
	    ipart._D = 0.0;                       /* set integer */
	    *odd = 0;
	 } else {
	    ipart._D = x;
	    if ( xexp <= L_FRAC ) {               /* integral bits are only in lower word */
	       ipart.st._LL = ( ipart.st._LL & ~( (unsigned long)MASK_FRAC >> xexp ));
	                                          /* clrear fractional bits */ 
	       ipart.st._LH = 0;
	       if ( (( ipart.st._LL & ((unsigned long)MASK_FRAC >> (xexp -1))) == 0 ) && ( xexp != 0 ) ) { 
	                                          /* 2^0 is LOW */
	          *odd = 1;
	       } else if ( xexp != 0 ) {
	          *odd = -1;
	       } else {                           /* return value is 1.0 */
	          *odd = -1;
	       }
	    } else {
	       ipart.st._LH = ( ipart.st._LH & ~( ~(unsigned long)0 >> ( xexp - L_FRAC )));
	                                          /* clrear fractional bits */ 
	       if ( (ipart.st._LH & (~(unsigned long)0 >> (xexp - L_FRAC -1))) == 0 ) {
	                                          /* 2^0 is LOW */
	          *odd = 1;
	       } else {
	          *odd = -1;
	       }
	    }
	 }
	 return( ipart._D );
	}


double
_Xfrac( x, pxexp )

	/************************************************************************/
	/*                                                                      */
	/*  The _Xfrac function devides x into the special number and an        */
	/*  integral power of 2.                                                */
	/*  The integer is stored in the long object pointed to by pxexp.       */
	/*  The argument x is positive number, not zero, Inf.                   */
	/*                                                                      */
	/*  RETURN    normal: the special number [ sqrt(2)/2 , sqrt(2) ]        */
	/*                      x = (return value) * 2^*pxexp                   */
	/*                                                                      */
	/************************************************************************/

/* ALGORITHM
 * 1.get the fraction part as xfrac and exponent part as *pxexp from the argument
 * 2.shift xfrac to the left to make hidden bit HIGH
 *   subtract shift number from *pxexp
 *     ( the argument is not normalized number )
 * 3.clear hidden bit   ( the argument is not normalized number )
 * 4.set exponent of xfrac
 * 5.if xfrac is less than sqrt(2)/2
 *      xfrac = 2 * xfrac
 *      *pxexp = *pxexp -1 
 */

	double         x;           /* positive and not zero,Inf number   */
	long           *pxexp;      /* the pointer, the integer is stored */
	{
	 DBL_LNG        xfrac;      /* the fraction part of x             */
	 const static double  h_sqrt2 = 0.70710678118654752440;

	 xfrac._D = x;                
	 *pxexp = (( xfrac.st._LL & MASK_EXP ) >> L_FRAC ); /* get exponent part  */
	                                          /* MASK_EXP = 0x7ff00000 */
	 xfrac.st._LL = ( xfrac.st._LL & MASK_FRAC ); /* get fraction part  */
	                                          /* MASK_FRAC = 0x000fffff */
	 if ( *pxexp == 0 ) {                     /* x is denormal number */
	    if ( xfrac.st._LL == 0 ) {            /* lower word is 0 */ 
	       xfrac.st._LL = xfrac.st._LH;       /* shift to the left by 1 word */
	       xfrac.st._LH = 0;
	       *pxexp = - ONE_WORD;
	    }
	    while ( xfrac.st._LL < DBL_MIN_L ) {  /* hidden bit is not HIGH */
	                                          /* DBL_MIN_L = 0x00100000 */
	       xfrac.st._LL = (( xfrac.st._LL << 1 ) | ( xfrac.st._LH >> ( ONE_WORD -1 ) ));
	       xfrac.st._LH = ( xfrac.st._LH << 1 ); /* shift to the left */
	       --*pxexp;
	    }
	    while ( xfrac.st._LL >= ( DBL_MIN_L << 1 ) ) { /* for overshoot */
	                                          /* shift to the right */
	       xfrac.st._LH = (( xfrac.st._LH >> 1 ) | ( xfrac.st._LL << ( ONE_WORD - 1 ) ));
	       xfrac.st._LL = ( xfrac.st._LL >> 1 );
	       ++*pxexp;
	    }
	    xfrac.st._LL = ( xfrac.st._LL & MASK_FRAC );  /* clear hidden bit */
	    ++*pxexp;                             /* x = demormal  -->  *pxexp <= 0 */
	 }

	 xfrac.st._LL = ( xfrac.st._LL | ( NUM1_L - DBL_MIN_L )); /* | 0x3feooooo */
	                                          /* set exponet */
	 *pxexp = ( *pxexp - ((unsigned long)( NUM1_L - DBL_MIN_L ) >> L_FRAC )); /* get off offset */
	 if ( xfrac._D < h_sqrt2 ) {              /* |xfrac| < sqrt(2)/2 */
	    xfrac._D = xfrac._D * 2.0;            /* sqrt(2)/2 <= xfrac < sqrt(2) */
	    --*pxexp;
	 }
	 return( xfrac._D );
	}



